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The International Mass Loading Service 

Leonid Petrov 

Abstract The International Mass Loading Service computes four loadings: 
a) atmospheric pressure loading; b) land water storage loading; c) oceanic 
tidal loading; and d) non-tidal oceanic loading. The service provides to users 
the mass loading time series in three forms: 1) pre-computed time series for a 
list of 849 space geodesy stations; 2) pre-computed time series on the global 
1° x 1° grid; and 3) on-demand Internet service for a list of stations and a 
time range specified by the user. The loading displacements are provided for 
the time period from 1979.01.01 through present, updated on an hourly basis, 
and have latencies 8-20 hours. 

1 Introduction 

Loading is a crustal deformation caused by a redistribution of air or water 
mass. In particular, loading caused by the redistribution of air mass is called 
atmospheric pressure loading, redistribution of continental water mass in a 
form of snow cover, soil moisture, and ground water causes continental water 
storage loading (sometimes also referred to as “hydrological loading”), redistri¬ 
bution of oceanic water mass causes ocean loading, which in turn is sub-divided 
into the ocean tidal and ocean non-tidal loadings. The loading deformation on 
average has the rms of 3 mm but may reach 60 mm (M 2 tidal ocean loading 
near British island). Calculation of mass loading is a computationally inten¬ 
sive procedure. For a case of ocean loading, the coefficients of displacement 
expansion can be computed once and forever. Computation of other loadings 
requires knowledge of surface pressure changes caused by a redistribution of air 
and water masses that is highly volatile and cannot be computed beforehand. 
Acquiring information about time series of global mass redistribution poses a 
serious logistical problem that impeded implementation of data reduction for 
mass loading. 

Realizing the magnitude of logistical problem, the mass loading service 
at NASA Goddard Space Flight Center was established in December 2002 
(Petrov and Boy 2004). The service was limited to the atmospheric pressure 
loading. It provided for the geodetic community the time series of 3D dis¬ 
placements of several hundred space geodesy stations, as well as the global 
displacement field at the 1° x 1° grid. The time series were updated daily 
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and had latency around 3 days. The service quickly became very popular and 
became the main source of information about atmospheric pressure loading. 

Recently, a decision was made to make a deep upgrade of the service. 
The upgrade targeted the following areas: a) to extend the service to loading 
caused by land water storage and non-tidal ocean loading; b) to support high 
resolution models of the atmosphere and land water storage; c) to take into 
account effects of local topography on surface pressure in mountain regions; 
d) to improve latency; e) to provide a user an ability to compute loading for 
user-selected stations on-demand. In the rest of the paper I will describe the 
approaches used for this upgrade. 


2 The use of high resolution models for loading computation 


The original atmospheric pressure loading service used the 2D NCEP Reanal¬ 
ysis surface pressure field (Kalnay et al. 1996) at a regular grid with a spatial 
resolution 2.5° x 2.5°. Modern models have much higher resolutions: for in¬ 
stance, the GEOS-FP model has resolution 0.3125° x 0.25°. The traditional 
approach for loading computation at a point with coordinate r involved a 
numerical evaluation of the integral of a convolution type (Farrell 1972): 
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where AP(r',t) is the pressure caused by mass redistribution, L((f>, A) — is 
the land-sea mask, the share of land in an elementary cell, and G(ip) are the 
Green’s functions defined as 
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The problem is that this algorithm has complexity G(d 1 2 3 4 ), where d is the 
spatial grid size, i.e. it grows very rapidly with an increase of spatial resolu¬ 
tion. It becomes impractical to use convolution for loading computation using 
models with a high spatial resolution. The alternative is to use the spherical 
harmonic transform approach. The algorithm involves the following steps: 


1. forming the pressure difference with respect to the average; 

2. transforming the surface pressure field to the regular grid with a higher 
resolution (upgridding): 2(D+1)+1 x 4(D+1) over latitude and longitude, 
where D is degree of the expansion; 

3. multiplying the surface pressure field with the land-sea mask defined as a 
share of land in a cell; 

4. spherical harmonic transform of degree/order D; 
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5. scaling the output of the spherical harmonic transform with Love numbers 
h' n and l' n of the corresponding degree n: 

3 h' 
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where p® is the mean Earth’s density and go is the equatorial gravity 
acceleration. The expression under the integral is the spherical harmonics 
() ™ of the the pressure field with the land-sea mask applied. 

6. inverse spherical harmonic transform: 
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This algorithm is equivalent to eqn (1) when D —^ oo, but it has complex¬ 
ity 0(d 3 ). It outperforms the convolution algorithm when D > 30. Numerical 
tests showed that in order to have errors in loading computation everywhere 
on the Earth less than 0.15 mm, degree/order 1023 is sufficient. It may sound 
counter-intuitive why such high resolution (0.088°) is needed, since the reso¬ 
lution of numerical models is one order of magnitude coarser. We should bear 
in mind that although the output of numerical models does not have signal 
at degree/order greater than 200-400, the product of the surface pressure and 
the land-sea mask is not band-limited and its spherical harmonic transform is 
not zero at any degree/order. 


3 Mass redistribution models 

Three numerical weather models developed at the NASA Global Modeling and 
Assimilation Office (GMAO) are used for loading computation: 

— MER.RA (Modern-Era Retrospective analysis for Research and Applica¬ 
tions) (Rienecker, et al. 2011). Resolution: 0.67° x0.5° x72 layersxG^, runs 
from 1979.01.01 through present, latency 20 d -60 d . This model is frozen and 
it is considered the most stable. 

— GEOS-FP (Global Earth Observing System Forward Processing) (Molod et al. 
2012). Resolution: 0.3125° x 0.25° x 72 layers x 3 /l , runs from 2011.09.01 
through present, latency 6 h -15 h . This is the operational model, updated 
approximately once a year. 
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— GEOS-FPIT (Global Earth Observing System Forward Processing Instru¬ 
mental Team) (Rienecker, et al. 2008). Resolution: 0.625°x0.5°x72 layersx 
3 h , runs from 2000.01.01 through present, latency 6 h -25 h . In terms in sta¬ 
bility this model is intermediate between MERRA and GEOS-FP, but it 
has a low latency. 

The surface pressure is computed from a 3D model. This process involves 
several steps. Firstly, each column of the output at the native, irregular , 
terrain-following grid is interpolated to the column at a new regular grid that 
is formally extrapolated down to -1000 m and up to 90,000 m. Then the at¬ 
mospheric pressure at a given epoch is expanded into the tensor product of 
B-splines over the entire Earth. Using the expansion coefficients, the pressure 
on the surface at resolution D1023 (0.088° x 0.088°) is computed. The height of 
the surface is derived from 30" x 30" GTOPO30 model 1 by averaging over cells 
of the D1023 grid. Using the expansion coefficients, the atmospheric pressure 
on that surface is computed. This procedure mitigates effects of orthography: 
in mountainous regions a node of a coarse grid may fall into a valley or a ridge 
and therefore, may not be representative for an average pressure of the cell. 
Three land water storage models are used for loading computation: 

— GLDAS NOAH025 (Global Land Data Assimilation System) (Rodell et al. 
2004). Resolution: 0.25° x 0.25° x 3A runs from 2000.01.24 through present, 
latency 35 d -70 d . 

— MERRA TWLAND (Rcichle et al. 2011). Resolution: 0.67°x0.5° x6 h , runs 
from 1979.01.01 through present, latency 35 d -60 d . This model is considered 
the most stable. 

— GEOS-FPIT TWLAND. Resolution: 0.625° x0.5°xl h , runs from 2000.01.01 
through present, latency 6 -25 . It was found that hourly time resolution is 
excessive for loading computation. The resolution was reduced to 3 hours. 

Upgridding involves refining the pressure field according to the fine land- 
sea mask. If a cell at the new 0.088° x 0.088° grid falls in the area that was 
ocean in the old grid, the pressure of the water equivalent of soil moisture 
and/or snow cover is computed by interpolation from surrounding cells that 
are land in the original grid with applying Gaussian smoothing. 

Non-tidal ocean loading is computed from the Ocean Model for Circulation 
and Tides (OMCT) (Thomas 2002). The original resolution of the model is 
1° x 1° x 6 h , latency: 10 d -60 d . However, the output of the original model is not 
available, only its spherical harmonic transform truncated at degree/order 100. 
Ugridding the OMCT model involves an iterative procedure that resembles the 
CLEAN algorithm used in radio astronomy for image restoration: it exploits 
the facts that the ocean bottom pressure is zero at land and the bottom 
pressure is relatively smooth in the ocean, except a jump to the zero at the 
shore. 

Two models of ocean tidal loading are used: the GOT4.8 (Ray 2013) and 
FES2012 (Carrere et al. 2012). They are upgridded to degree/order 2047 in a 
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similar way as it was done for land water storage, except reversal of land and 
sea cells. 


4 Processing pipeline 


The two servers of the international mass loading service that work indepen¬ 
dently check every hour whether new data appeared. If the new data appeared, 
they are downloaded, decoded, up-gridded, and the surface pressure anomaly 
at the D1023 grid is computed by subtracting a model that includes the mean 
surface pressure value, sine and cosine amplitudes of pressure variations in a 
range of frequencies in the diurnal, semi-diurnal, ter-diurnal and four-diurnal 
bands. Then the spherical harmonic transform of degree/order 1023 of the 
pressure field anomaly is computed and scaled by Love numbers of the corre¬ 
sponding order. The coefficients V™ and H™ in eqn (3) are stored. They are 
used for loading computations in three ways: 


1. Computing loading at the D89 grid (1° x 1°). This is done in the following 
way: the spherical harmonic transform of degree/order D1023 is padded 
with zeroes to degree/order D1079. The coefficients V™, H™ are underwent 
the inverse spherical harmonic transform and produce the loading field in 
local Up, East, and North direction at the D1079 grid (1/12° x 1/12°). 
Every 12th element of the intermediate D1079 grid is written in the output 
file. 

2. Computing loading for a set of 849 commonly used GNSS, SLR, DORIS, 
and VLBI stations. 

3. Computing loading on-demand for the set of stations supplied by the user. 
A user fills the Web form where he or she specifies the model, the range of 
dates and the list of stations with their Cartesian coordinates. When the 
loading computation is finished, a user can retrieve the files with results. 


The loading displacements are computed using the Love numbers defined 
in the coordinate system with the origin at the center of mass of the total 
Earth: the solid Earth and the fluid under consideration. For some applica¬ 
tions displacements with respect to the center of mass of the solid Earth are 
desirable. The loading Love numbers differ between these systems only for de¬ 
gree 1. The International Mass Loading Service computes the differential load¬ 
ing displacements between these two systems. Such a displacement, called the 
“degree one displacement” uses only two terms of degree 1 in expression (4): 

1 + h i ,/m „„ A um _ 1 + l l 
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—j-f — VT and H[ n = —n- 
n 1 L l 

is added to the displacement with respect to the center of the total mass, the 

sum is the displacement with respect to the center of mass of the solid Earth. 
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5 Validation 

VLBI observations for the period of 2001.01.01 - 2014.07.01 were used for 
loading validation. The same technique was applied as we used for loading 
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Table 1 Estimates of admittance factors for Up (UP), East (EA) and North (NO) compo¬ 
nents for three different loading models from the global least squares solution using geodetic 
VLBI group delays. 


Atm GEOS-FPIT UP 

0.963 

± 

0.023 

Atm GEOS-FPIT EA 

0.609 

± 

0.049 

Atm GEOS-FPIT NO 

1.027 

± 

0.041 

Lws GEOS-FPIT UP 

0.955 

± 

0.016 

Lws GEOS-FPIT EA 

0.804 

± 

0.029 

Lws GEOS-FPIT NO 

0.886 

± 

0.024 

Lws NOAH025 UP 

1.220 

± 

0.013 

Lws NOAH025 EA 

0.660 

± 

0.030 

Lws NOAH025 NO 

0.826 

± 

0.033 


validation in Petrov and Boy (2004): the global admittance factors were es¬ 
timated from the data together with estimation of site positions, velocities, 
the Earth orientation parameters, source coordinates and nuisance parame¬ 
ters such as clock functions and atmosphere path delays in zenith direction 
(see Table 1). The partial derivative for admittance factors was the contribu¬ 
tion of the loading displacement into path delay. If the model is perfect, the 
admittance factor will approach to unity. 

Surprisingly, the GEOS-FPIT land water storage model has admittance 
factors closer to unity than the GLDAS NOAH025 model. This is important for 
practical applications, since the GEOSFPIT model has much lower latencies 
than the GLDAS NOAH025 model. 


6 Using the International Mass Loading Service 

The gridded loading displacements are useful for visualization of the loading 
field and for computation of integrals over the area. However, a user should 
be aware that the held of loading displacement near the coastal area is not 
smooth. Therefore, using gridded loading for data reduction by interpolation 
the displacement held to the position of a given station may cause significant 
errors. This problem is illustrated in Figures 1-2 for a case of ocean loading 
near Newfoundland. The M 2 ocean loading displacement has the vertical am¬ 
plitude ~ 30 mm, but interpolation errors exceed 30% within 100 km of the 
coastal area when the 1.0° x 1.0° grid is used. The errors are in excess of 30% 
within 30 km from the coast when the 0.25° x 0.25° grid is used. They fall 
below 1 mm only when the grid with a resolution 0.05° x 0.05° or hner is used. 

Gridded loading at 1° x 1° or 0.25° x 0.25° resolutions should never be 
used for data reduction. The International Mass Loading Service computes 
loadings for 849 stations directly without the use of interpolation. This is suf¬ 
ficient for processing SLR, DORIS, and VLBI observations, since new stations 
are introduced infrequently. New GNSS stations are introduced much more 
frequently, and the situation when the GNSS station of interest is absent from 
the list of stations with pre-computed loading is more common. Figure 3 shows 
the Web interface that implements the on-demand loading computation. The 
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Fig. 1 Mass loading caused by the M 2 ocean tide near Newfoundland island computed 
with two resolutions: 0.01° X 0.01° grid (Left) and 1.0° X 1.0° grid (Right) 



Fig. 2 The difference of mass loading caused by the M 2 ocean tide computed with two 
resolutions: 1.0° X 1.0° grid versus 0.01° X 0.01° grid (Left) and 0.25° X 0.25° grid versus 
0.01° X 0.01° grid (Left) and 



on-demand computational procedure uses V™'(t) and H™(t) coefficients that 
are evaluated and stored as soon as new data arrive. For this reason, the on- 
demand procedure is relatively quick: computation of loading displacements 
for 20 stations for a 1 year interval with the time step of 3 hours, in total 
58,400 loading displacements, takes 20 minutes. 

7 Conclusions and future work 

At present, the International Mass Loading Service offers to the geodetic com¬ 
munity computation of 3D displacements caused by the atmospheric pres¬ 
sure loading, land water storage loading, tidal and non-tidal ocean loading, 
free of charge, 24/7 with a latency from 8 hours (atmospheric and land wa¬ 
ter storage loading) to 30 days (non-tidal loading). The URL of the pri¬ 
mary server is http://massloading.net, the URL of the secondary server 
is http://alt.massloading.net. The loading displacement were validated 
against the dataset of global VLBI observations for 2001 2014. 

Further development: a) using weather forecast to 0-24^ in the future. 
Latency will be eliminated. Accuracy degradation with respect to an assimila¬ 
tion model: 20% for the current instant; b) using OPeNDAP protocol for data 
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Fig. 3 User interface to the on-demand computation of mass loading displacements 

<- -* o *- |oTurbo massloading.net/atm/_★ [ a - Search with GoT 

Compute displacements caused by atmospheric pressure loading on-demand * 

You can order computation of the time series for the stations of your interest. You need to prepare a station file in plain ascii 
that has four columns separated by one or more blanks: 

Station_name X- coordinate Y-coordinate Z-coordinate 

Station name should have no more than 8 characters. X,Y,Z are Cartesian coordinates of the station of interest in a crust-fixed 
coordinates system. Units are meters. Here is an example . 

O MERRA from 19790101_0000 through 20140731_1800 
O GEOS507 from 20110901_0000 through 20130614 0300 
Model: O GEOS511 from 20130611_0000 through 20140820_0300 

O GEOSFP from 20140801_0000 through 20140929_0900 
* GEOSFP1T from 20000101_0000 through 20140929_0900 
Mode: * Time series O harmonics coefficients 

Frame: # Center of mass o Center of figure 

Start date: |20i4.06.06 | Format: YYYY.MM.DD_hh:mm:ss 

Stop date: |20i4.07.n ~] Format: YYYY.MM.DD_hh:mm:ss = 

Station file: I'ytmpQ.lnp' IfcFioose... I 

E-mail address: I I Optional 

| Submit") 

Results of on-demand computations are accessible from here . 

Plots of displacements caused by harmonic variations of atmospheric pressure 

distribution; c) automation of loading displacement ingestion: development 
of a client library that communicates with the loading servers automatically; 
d) generation of time series of the sum of all loadings on the fly with a user- 
requested time step; e) computing loading due to water level changes in lakes 
and big rivers. 

This project was supported by NASA Earth Surface and Interior program, 
grant NNX12AQ29G. 
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